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ABSTRACT 

Models of the Sagittarius Stream have consistently found that the Milky Way disc is 
O . oriented such that its symmetry axis is along the intermediate axis of the triaxial dark 

matter halo. We attempt to build models of disc galaxies in such an "intermediate-axis 
^ , orientation" . We do this with three models. In the first two cases we simply rigidly grow 

a disc in a triaxial halo such that the disc ends up perpendicular to the intermediate 
axis. We also attempt to coax a disc to form in an intermediate-axis orientation by 
producing a gas-|-dark matter triaxial system with gas angular momentum about the 
^ ' intermediate axis. In all cases we fail to produce systems which remain with stellar 

angular momentum aligned with the halo's intermediate axis. For one of these unstable 
simulations we show that the potential is even rounder then the models of the Milky 
Way potential in the region probed by the Sagittarius Stream. We conclude that the 
Milky Way's disc is very unlikely to be in an intermediate axis orientation. However 
we find that a disc can persist off' one of the principal planes of the potential. We 
' propose that the disc of the Milky Way must be tilted relative to the principal axes of 

, the dark matter halo. Direct confirmation of this prediction would constitute a critical 

test of MONO. 

Key words: Galaxy: evolution - Galaxy: formation - Galaxy: halo - Galaxy: kine- 
' matics and dynamics - Galaxy: structure - galaxies: haloes 
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1 INTRODUCTION 

Dark matter hal os in purely coll i sionle s s simulations are gen- 
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in both the type and shape of orbits (|Valluri et al.ll20ig ) 
Nonetheless halos remain triaxial at large radii. 

The Sagittarius dwarf tidal stream, which extends to 
~ 60 kpc from the Galactic centre, has been used to con- 
strain the shape of the Milky Way's halo with varying re- 
sults. Noting that the tidal debris is distributed on a great 
circle, llbata et al.l fcoOl) concluded that the halo is nearly 
spherical. Likewise IPellhauer et al.l (|2006l ) argued that the 
position of the bifurcation in the tidal stream, which they in- 
terpreted as two wraps o f the stream, can be explained if the 
halo is close to spherical. iMarti'nez-Delgado et al.l (|2004l ) and 



2 Debattista et al. 



Ijohnston et al.l (|2005l ) instead found a mildly oblat e halo 

i c/g ^ 0.9) flattened in the same sense as the disc. iHelmil 
2004ah meanwhile argued that the trailing part of the Sagit- 
tarius stream is too dynamically young to pro vide a s t ringent 
constraint. Using instead the leading stream, iHelmil l|2004bl ) 
found evidence for a prolate halo with c /a ^ 0.6 and with 
its long axis perpendicular to the disc. iLaw et al] (|2009l ) 
were the first to demonstrate that simultaneously fitting the 
density and kinematics of the Sagittarius stream requires a 
triaxial (rather than oblate/prolate) potential. A surpris- 
ing property of this potential is that its intermediate axis 
is aligned with the short axis of the disc (a relative orienta- 
tion we will refer t o as t he "intermediate-axis orientation"). 
iLaw fc MaiewskH (|2010l . hereafter LMIO) presented a suite 
of ~ 500 A'^-body simulations of the tidal disruption of the 
Sagittarius dwarf in a fixed potential. The simulations were 
compared with a large number of constraints including (i) 
the position and velocity of the Sagittarius dwarf with its 
velocity vector in the orbital plane of the young trailing tail, 
(ii) the radial velocity and velocity dispersion in the trail- 
ing tidal tail and (iii) the angular location, width and radial 
velocities of the leading tail. The best fitting model is in 
the intermediate-axis orientation with (6/a)a> = 0.99 and 
{c/a)g, = 0.72 between 20 kpc and 60 kpc, and with the 
m ajor and minor axes in the plane of the disc. An analysis 
bv lDeg fc Widrowl l|2013 !) that varies also the parameters of 
the bulge-|-disc of the Milky Way still finds a disc in the 
intermediate-axis orientation. Recent extended mapping of 
the Sagittarius stream in the southern Galactic hemisphere 
finds a stream consistent with the LMIO model (|Slater et al.l 
I2OI2I ). 

Triaxial potentials are populated by box orbits (which 
get arbitrarily close to the centre of the potential) and 
tube orbits (which have a fixed sense of rotation rela- 
tive to one of the principal axes). The stability of tube 
orbits about each of the principal axes of a tria:xial po- 
tential has been studied extensively: tube orbits are sta- 
ble around the short and long axes, but not around 
the intermediate axis (Heiligman & Schwarzschild' '1979"; 
I Goodman & Schwarzschild 1981: Wilkinson & James 1982), 



even when planar ( Adams et al.ll2007[^ Carpintero fc Muzziol 
|2012| ). In semi-cosmological models. lAumer fc White! (|2013l ) 
showed that discs are most stable when the angular momen- 
tum is aligned wit h the minor axis of the halo. The model of 
iLaw fc Maiewslq (i2010) therefore challenges the view that 
the instability of intermediate-axis tube (lAT) orbits pro- 
hibits discs from forming in this orientation. One way in 
which this discrepancy might be resolved is if in the vicinity 
of the disc it dominates the net potential, which becomes 
fiattened like the disc. Then the near-circular orbits in the 
disc are orbiting around the local short axis, and therefore 
in a stable configuration (K. Johnston, private communica- 
tion) . 

In this paper we show that discs are unable to persist 
in an intermediate-axis orientation. We use both simulations 
in which discs are grown inside isolated triaxial halos as well 
as a simulation of a galaxy forming out of gas with angular 
momentum about the intermediate axis of a triaxial halo. 
In Section [2] we discuss the methods used in this paper, 
including the initial conditions of the stars, dark matter and 
gas. Section[3]presents the evolution of the models. We draw 
our conclusions in Section [l] The Appendix presents our 



interpretation for why the intermediate-axis orientation is 
unstable based on an orbital study. 



2 NUMERICAL METHODS 

2.1 Constructing collisionless initial conditions 

As in iDebattista et al.l (|2008l '). we formed triaxial halos via 
the merger of three or more spherical halos (^ Moore et al.l 
i2004i) . The mergers, and all subse quent collisio nless simula- 
tions, were evolved with pkdgrav (|Stadel200ll '). an efficient, 
multi-stepping, parallel treecode. The spherical halos were 
generated from a d istrib ution function using the method of 
iKazantzidis et al 

I 12004") with each halo composed of two 
mass species arranged on shells. The outer shell has more 
massive particles than the inner one, inc reasing the effec- 
tivc resolution in the centre. As shown in IDebattista et al.l 
(2008), a large part of the particle mass segregation persists 
after the mergers and the inner region remains dominated 
by the lower mass particles. 

We produced two dark-matter-only triaxial halos, which 
wc refer t o as A and C; halo A was presented already in 
Debattist a et al.l (|2008l ). These halos were constructed from 
two mergers. In both cases the first merger placed two con- 
centration c = 10 halos 800 kpc apart approaching each 
other at 50 kms^^, producing a prolate merged halo. Halo 
A was generated by the head-on merger of two copies of this 
halo starting at rest 400 kpc apart. For halo C, after the 
first merger, a third halo, with c = 20, was merged from 
100 kpc along the first remnant's minor axis. The top two 
panels of Figure [T] plot the shap e and trictxiality of thes e two 
halos, measured as described in lDebattista et aL 1(2008^, be- 
fore any discs are introduced. The triajciality param eter is 
defined as T = (a^ ~ b'^)/{a'^ - c^) (|Franx et al.lll99ll '). Halo 
A is highly prolate but has only a mild triaxiality T ~ 0.9; 
its shape however is very constant out to 100 kpc. Halo C 
instead has a radially varying T ranging from ~ 0.9 at the 
centre to ~ 0.3 at 50 kpc. Halo C is considerably rounder 
than halo A everywhere within the inner 100 kpc. Table [1] 
lists the properties of the haloj^. 

The outer particles are ~ 19x more massive than the 
inner ones in halo A. Halo C has four particle mass species, 
with the outer particles up to 19 x more massive than the 
inner ones. Halo A has 4 million particles while halo C has 3 
million. We used a softening parameter e = 0.1 kpc for low 
mass particles and e — 0.5 kpc for the more massive species. 

Once we produced the triaxial halos, we inserted a disc 
of particles which initially remained rigid. The disc distri- 
bution was, in both cases, exponential with scale-length 
Rd = 3 kpc and Gaussian scale-height ~ 0.05-Rd. The 
discs are comprised of 300K equal-mass particles each with 
a softening e = 60 — 100 pc. Initially the disc has negligible 
mass but this grows linearly over 5 Gyr. During this time, 
the halo particles are free to move and achieve equilibrium 
with the growing disc. 

The disc in halo A is grown to a mass of 1.75 x 10^^ M0. 



^ We use a differe nt c onvention from IDebattista et "aL I l l2008f) . 
IValluri et al] 1 I2OIOI) and lValluri et al.l 1I2OI2 I'). who used the radius 
at which p = 200pcrit- Here r2oo is the radius within which the 
enclosed mass has average density 200pcrit- 
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Table 1. The halos used in the simulations. The properties listed 
are for the halo after the last merger and before the discs have 
been grown. Np and Ng are the number of dark matter and gas 
particles within r2oo i and ^200 is the halo mass within the virial 
radius, r200. Density axes-ratios b/a and c/a are by-eye averaged 
over the inner 20 kpc (see Figure [TJ. 



The disc is placed in an intermediate-ajcis or ientation and 
we therefore refer to this model as model lAl. IValluri et al.l 
(|2012l ) presented an orbital analysis of the halo in this model 
at t = 0; there the model is also referred to as lAl. For some 
of our analysis, we also present a version of this model with 
the disc at a mass of only 7 x 10^° Mq, which we refer to as 
model IA2. This evolves quite similar to model lAl. The disc 
in halo C instead is placed with its symmetry axis along the 
halo's long axis, so we refer to it as model LCI. This disc has 
a final mass of 1.4 x 10^^^ Mq . The high disc masses in models 
lAl and LCI allow us argue that even a high mass does not 
offer a disc protection against the unstable intermediate- 
axis orientation, but we have checked that lower mass discs 
(including in IA2) are also unstable in this orientation. 

We set the kinematics of the final discs to give 
constant zg and Toomre-Q = 1.5, as described in 
iDebattista fc SellwoodI (|2000! V For this we calculate the po- 
tential using a hybrid polar-grid code with the disc on a 
cylin drical grid and the halo on a spherical grid (jSellwoodl 
l2003l ). In setting the disc kinematics, we azimuthally average 
radial and vertical forces; thus our discs are initially not in 
perfect equilibrium. Equilibrium is quickly established once 
the disc particles are free to move. In these simulations t = 
corresponds to the time at which we set the disc kinemat- 
ics. PKDGRAV is a multi-stepping tree code, with timesteps 
refined such that St = At/2" < ri{t/aY^'^, where e is the 
softening and a is the acceleration at a particle's current 
position. We use base timestep At = 5 Myr, rj — 0.2 and set 
the opening angle of the treecode to 6 = 0.7 in all cases. 



2.2 Initial conditions with gas 

We also present a simulation of a disc forming out of gas ro- 
tating about the intermediate axis of a triaxial halo , whic h 
we refer to as model Gil. As did lAumer fc White! (|2013f ). 
in our initial experiments we found that arbitrarily insert- 
ing rotating gas halos within pre-existing non-spherical dark 
matter halos leads to a substantial loss of gas angular mo- 
mentum. Our approach therefore is to include the gas, which 
is not allowed to cool, right from the start while merging ha- 
los to produce the triaxial system. We first set up a prolate 
halo with an equilibrium gas distribution by merging two 
spherical NFW dark matter halos as before. Each of the 
spherical initial halos has an embedded spherical hot gas 
component containing 10% of the total mass and following 
the same d ensity distribution. The initial halos have been 
described in lRoskar et al.l l|2008l ): each dark matter halo has 
a mass within the virial radius of 10^^ Mq. A temperature 
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Figure 1. Shape of dark matter halos A (top), C (middle) and 
Gil (bottom) before any of the discs/star formation are intro- 
duced. Solid, dashed and dotted lines show b/a, c/a and T, re- 
spectively. 



gradient in each halo ensures an initial gas pressure equilib- 
rium for an adiabatic equation of state. Gas vel ocities are ini- 
tialised to give a spin par ameter of A = 0.039 (jBuUock et al.l 
I2001I : Im accio etaP 120071 ). with specific angular momentum 
j oc R, where R is the cylindrical radius. Each halo used 10® 
particles in each of the gas and dark components. Gas parti- 
cles initially have masses 1.4 x 10^ Mq and softening 50 pc, 
the latter inherited by the star particles, while dark matter 
particles come in two mass flavors (10® Mq and 3.5 x 10® Mq 
inside and outside 200 kpc respectively) and with a soften- 
ing of 100 pc. The two halos are placed 500 kpc apart along 
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Gil 



Figure 2. Cartoon representation of the merger geometry that 
produces a triaxial halo with gas angular momentum along the 
intermediate axis, model Gil. The (red) filled arrows indicate 
the orientation of the gas angular momentum relative to the two 
merging prolate halos while the (black) open arrows indicate the 
relative motion of the 2 halos. 

the i-axis and are initially moving towards each other at a 
relative velocity of 100 kms^^. 

After the first merger the resulting halo is prolate, elon- 
gated along the i-axis, with (c/o) ~ 0.65 and angular mo- 
mentum along the short (z) axis. We produce a triaxial halo 
by merging two copies of this prolate system (for a total of 
4 X 10*^ particles in each of the gas and dark matter com- 
ponents). In order to align the gas angular momentum with 
the intermediate axis of the halo we first rotate the prolate 
system about the long axis so the angular momentum vec- 
tor is along the y-axis, then rotate two copies of the prolate 
halo about the z-axis by -1-30° in one case and by —30° in 
the other. This merger geometry for the two prolate halos 
is illustrated in Figure [S] Merging these two halos from a 
separation of 500 kpc along the a;-axis with a relative veloc- 
ity of 100 kms~^ produces a quite prolate halo with only a 
very mild triaxiality (T ^ 0.93 within the inner 100 kpc), as 
shown in the bottom panel of Figure [T] 

This sim ulation was evolved with GASOLINE 
jWadslev et al.l [2004 ). the smooth particle hydrody- 
namics (SPH) version of pkdgrav. We use a base timestep 
of 10 Myr with a refinement parameter r] = 0.175. During 
the mergers, and for some time after, we evolve the sys- 
tem adiabatically without gas cooling or star formation. 
After, we switch on gas cooling, star formation and stella r 
feedback using the prescriptions of IStinson et al.l l|2006l ). 
A gas particle undergoes star formation if it has number 
density n > O.lcm"^, temperature T < 15,000 K and is 
part of a converging fiow; efficiency of star formation is 
0.05, i.e. 5% of gas particles eligible to form stars do so 
per dynamical time. Star particles form with an initial 
mass of 1/3 that of the parent gas particle, which at our 
resolution corresponds to 4.6 x 10* Mq. Gas particles can 
spawn multiple star particles but once they drop below 
1/5 of their initial mass the remaining mass is distributed 
amongst the nearest neighbours, leading to a decreasing 
number of gas particles. Each star particle represents an 
entir e stellar population with a Miller-Scalo (j Miller fc S caloj 
1X9791 ) initial mass function. The evolution of star particles 
includes AGB stellar winds and feedback from type II and 
type la supernovae, with their energy injected into the 
interstellar medium (ISM) . Each supernova releases 4 x 10^" 
ergs into the ISM. The effect of the supernovae explosions is 
modelled at the su b-grid level as a bla st-wave propagating 
through the ISM (|Stinson et al.ll2006l ). We again use an 
opening angle oi 9 = 0.7. The timestep of gas particles also 
satisfies the condition Stgas = r)couranth/[{l + a)c + Pfimax], 



where rjcourant = 0.4, h is the SPH smoothing length, a is 
the shear coefficient, which is set to 1, /3 = 2 is the v i scosit y 
coefficient and fimax is described in IWadslev et al.l (|2004h . 
The SPH kernel is defined using the 32 nearest neighbours. 
Gas cooling is calculated without taking into account the 
gas metallicity. These prescriptions have been shown to lead 
to re alistic Milky Way-type galaxies (jRoskar et al.l |2012| . 
l2013h . In this run, t = corresponds to the time at which 
gas cooling is switched on and star formation commences. 



2.3 Briggs figures 

We us e Briggs figu res, originally introduced for studying 
warps (|Briggsl [l990l ). to illustrate disc tilting in the simu- 
lations. A Briggs figure is a 2-D polar coordinate represen- 
tation of the direction of vectors. We decompose the stellar 
discs into 5 concentric rings of equal width extending to a 
radius of 15 kpc and for each ring plot the direction of the 
angular momentum vector in 2-D cylindrical polar coordi- 
nates. The tilt of the angular momentum vector from some 
fiducial 2-axis, 6, is plotted as the radial coordinate while the 
angle from some fiducial a;-axis, (j), is plotted as the angle co- 
ordinate. Briggs figures are useful for showing the evolution 
of disc orientation provided that the axes with respect to 
which the angles 9 and (p are defined are kept fixed. Note 
that the Briggs figure of a uniformly tilting disc consists of a 
set of coincident points, indicating that the angular momen- 
tum of the disc is everywhere aligned. A differentially tilting 
(i.e. warped) disc instead is represented by non-coincident 
points. In the coUisionless simulations we always set the z- 
axis to be the direction of the angular momentum of the 
initial disc. Note that th i s is d ifferent from the convention 
adopted bv lValluri et al.l (|2012l ). 



3 RESULTS 

3.1 Models lAl and IA2 

In models lAl and IA2 the i-axis is the pre-disc halo long 
axis while the y-axis is the short axis. Once the disc is grown, 
however, the inner halo is fiattened so the disc's vertical {i.e. 
symmetry) axis becomes the shortest axis of the inner halo. 
At larger radii the x and y axes continue to be the long 
and short axes of the halo, so we use these to specify the 
axes ordering. We compute the potential in the x = and 
y — planes, from which we measure the axes-ratios of the 
potential. The top panel of Figure [4] plots the equipotential 
axis ratios x^/z<i, and y^/z<i, for lAl. The pre-disc ordering 
of the axes of the potential is /z<i, > 1 > y^ j z^ ^ > 
Zij, > y<s> but after the disc is grown, within 20 kpc this 
becomes a;* > yi> > z^. 

The evolution of run lAl is shown in FigureO The disc 
tilts by 90° out of the initial plane within 4 Gyr. During this 
rapid tilting phase the disc does not warp substantially or 
precess (which can be seen from the fact that disc symme- 
try axis does not circulate about any axis). At f = 4 Gyr the 
disc has not yet settled, having overshot the minor axis ori- 
entation to ~ 120°. After 4 Gyr the disc precesses about 
the short axis while slowly settling into a short-axis orien- 
tation. Throughout this evolution, total angular momentum 
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Figure 3. Briggs figures (see Section 12.31 for a explanation of 
these figures) showing the evolution of run lAl at 1 Gyr inter- 
vals. Dotted circles are spaced at 20° intervals, with the outer 
solid circle corresponding to 6 = 120°. The centre of the disc is 
indicated by the (red) open circle while the remaining disc annuli 
are indicated by (red) crosses. The open (blue) star, square and 
triangle symbols indicate the direction of the pre-disc halo long, 
intermediate and short axes, respectively. 

is conserved to better than 1.5%, with angular momentum 
exchanged between the disc and the halo. 

Intermediate-axis tub e (lAT) orbits are unstable (e.g. 
iBinnev fc Tremainell200^ . p. 263). If the disc is perpendic- 
ular to the intermediate axis of the potential, then its stars 
would be on lAT orbits, which would render them unstable. 
After the disc has grown the net potential becomes so verti- 
cally flattened that the z-axis becomes the shortest axis of 
the potential in the disc's vicinity. This is the case also if just 
the halo potential is considered. Therefore the disc tube or- 
bits are stable because they are cocooned inside a vertically 
flattened halo and circulate about the shortest axes of their 
local potential. We confirm this by repeating the simulation 
with the halo particles frozen in place in model IA2. Then 
the disc does not tilt during 5 Gyr. 

The instability must therefore reside in the halo. In the 
Appendix we present evidence that the instability is driven 
by the response of tube orbits to a potential with a radially 
varying orientation. Because the halo has negligible angular 
momentum, it tilts without precessing, shepherding the disc 
along with it. Evidence that the halo is driving the tilting 
of the disc comes also from the small angular displacement 
between the disc and the inner halo. Close examination of 



in 




R [kpc] 

Figure 4. Equipotential axis ratios in models lAl (top) and LCI 
(bottom). The solid lines show x^/z^ while the dashed lines show 
y^/z^. The thick blue lines correspond to the halo before the disc 
is grown. The black and red lines show the full and halo potential 
shape after the disc is grown. The dotted horizontal lines indicate 
an axis ratio of unity. The z-axis is perpendicular to the disc at 
t = 0. 

Figure [3] shows that during the tilting phase (2-4 Gyr), for 
the disc is not the same as that for the halo minor axis. In 
Figure O the red points mark the direction of the disc angu- 
lar momentum; thus the disc orientation during the tilting 
phase is ahead of (larger <j!>) the great circle between the 
intermediate and short axes, along which the halo tilts. In 
order to demonstrate this, we again use the lower disc mass 
model IA2, since this distorts the inner halo to a lesser ex- 
tent. Other than the disc tilting more rapidly, this lower 
mass run evolves similar to run lAl. Figure [5] shows the 
evolution of the direction of the inner halo (solid lines) and 
of the disc (dashed lines) minor axes separately, by plotting 
the tilt angle 6 from the z-axis and position angle, (j) from the 
a;-axis. The disc and halo tilt away from the original vertical 
axis together, but the halo (j> is clearly closer to (j) — 90°, 
corresponding to the outer halo minor axis, than is the disc 
(j). Since the halo tilts almost directly into the minor axis, 
the disc angle (j> relative to the minor axis can be understood 
as the disc tilt relative to the halo needed to generate the 
gravitational torque needed to tilt the disc. Once the inner 
halo has settled, the misalignment between the disc and the 
inner halo leads to the damped precession seen after 4 Gyr. 
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Figure 5. Evolution of the relative orientation of the disc and 
inner halo in model IA2. The disc tilts rapidly from its initial 
orientation. The solid and dashed lines show the halo and disc 
orientations, respectively, at different times, as indicated at bot- 
tom. The black, blue, green, cyan and red lines indicate t = 0, 
1, 2, 3, and 3.5 Gyr, respectively. The z axis relative to which 9 
is measured is perpendicular to the initial disc, while the x-axis, 
which defines = 0, is the long axis. 



Since the instability is due to the halo, no matter how mas- 
sive the disc becomes (the halo-to-disc mass ratio within 15 
kpc is 1.6), this orientation can never be stable. 



3.2 Model LCI 

Before the disc is grown in run LCI, the direction vertical to 
the disc is the long axis of halo C. Figure |4] shows that the 
ordering of the axes is xq, at this stage, but once 

the disc is grown, the halo at r < 10 kpc flips by 90° , so that 
the intermediate axis becomes the axis orthogonal to the 
disc. The combination of the disc and halo potential then 
has j/$ > > z$ within 15 kpc, but beyond that point 
> Zcj > xj.. Thus while most of the disc is immersed 
perpendicular to the short axis of the local potential, at 
larger radii the disc symmetry axis is along the intermediate 
axis of the potential. Although the flip in the principal axes 




Figure 6. Briggs figures showing the evolution of run LCI at 
1 Gyr intervals. Dotted circles are spaced at 20° intervals, with 
the outer solid circle corresponding to 9 = 120°. The centre of 
the disc is indicated by the (red) open circle while the remaining 
disc annuli are indicated by (red) crosses. The open (blue) star, 
square and triangle symbols indicate the direction of the pre-disc 
halo long, intermediate and short axes, respectively. In the inner 
halo y > z > X while in the outer halo z > y > x once the disc is 
grown. 



of the density extends only to the inner halo, the flip in the 
axes of the potential extends till at least 80 kpc. The halo flip 
is probably related to the accretion history of halo C which 
included an accretion along the minor axis of a prolate halo. 
Indeed the inner halo major axis flips into the direction of 
the original accretion event. 

Thus in run LCI the outer part of the disc starts out 
immersed in a potential which, locally, has its intermedi- 
ate axis perpendicular to the disc. The disc in run LCI 
tilts very rapidly, initially towards the original intermediate- 
axis orientation and then dropping into a nearly short-axis 
orientation, as shown in Figure (6] The tilting rate reaches 
~ 30° Gyr~^ between 2 and 4 Gyr. This rapid, direct tilting 
is not accompanied by precession or warping. When we re- 
run the simulation with the halo frozen, the outer disc still 
tilts and forms a polar ring, while the inner disc tilts by only 
~ 15°. Thus lAT orbits of stars in the outer disc region are 
highly unstable. However the entire disc is not tilting be- 
cause of this instability. Given the lack of precession when 
the disc is live, we conclude that the inner halo of run LCI 
is also in an unstable orientation, much as in run lAl. 

Figure [7] shows the shape of the potential out to 80 
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Figure 8. Briggs figure for the gas within the inner 100 kpc in 
model Gil at i = 0, before any star formation. Dotted circles are 
spaced at 20° intervals, with the outer solid circle correspond- 
ing to S = 100°. The centre of the gas halo is indicated by the 
(red) open circle while the remaining shells are indicated by (red) 
crosses. Each shell is 10 kpc wide. The open (blue) star, square 
and triangle symbols indicate the direction of the pre-disc halo 
long, intermediate and short axes, respectively. The (black) filled 
star represents the orientation of the total gas angular momentum 
within this volume. 



kpc. The middle panel shows the radial profile of the po- 
tential axis ratios, 2;#/y$ and z^/j/jj. The longest axis of 
the potential is the y-axis (recall that the axis vertical to 
the disc is z). Beyond ~ 15 kpc, the potential intermedi- 
ate axis is the z-axis (i.e. perpendicular to the initial disc) 
and its shape, while not constant, does not vary substan- 
tially with radius. The bottom panel of Figure [7] compares 
the shape of mod el LCI and the Milky Way potent ial in the 
iLaw et all (|2009l ). LMIO and lDeg fc Widrowl (|2013l ) models. 
Model LCI has larger (i.e. rounder) (c/a)$ than all these 
models, while (b/a)^ is comparable to the best LMIO and 
iDeg fc Widrowl (|2013| ) TP models. For a spherical potential, 
(6/a)| + {c/a)% = 2; we measure deviation from sphericity 
as ^ = 2 — {b/a)% — (c/a)|. The bottom panel of Figure [7] 
plots contours of ^ which clearly shows that the potential 
in LCI is more nearly spherical in this region than are the 
Milky Way models. The instability of model LCI is therefore 
very probably shared by all these Milky Way models. 



(b/a). 

Figure 7. Top: Profile of axes ratios of the potential at t = in 
models lAl (thick lines) and IA2 (thin lines). The solid (dashed) 
line shows y^/x^ (z^/x^). Middle: Profile of axes ratios of the 
potential at t = in model LCI. The solid (dashed) line shows 
^s/y* (2* /;/#)■ Bottom: potential axes ratios of LCI (black cir- 
cles), lAl (blue circles) and IA2 (green circles) in the radial range 
16 ^ r/kpc ^ 60 (with the open blue circle showing 16 kpc and 
the filled circles showing larger radii) and Milky Way models. 
Dashed lines are contours of deviations from sphericity as de- 
fined in the text. 



3.3 Model Gil 

Figure [8] shows the initial angular momentum of the gas 
within the inner 100 kpc of model Gil. The total angular 
momentum within this region is very well aligned with the 
intermediate axis of the halo. Only within 30 kpc is the gas 
angular momentum not in this orientation, but this corre- 
sponds to a tiny fraction of the total angular momentum of 
this gas. 

During the first 2 Gyr of evolution after gas cooling and 
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Figure 9. The stellar disc in run Gil at 6 Gyr as seen in its 
intrinsic planes: face-on (bottom-left) and two orthogonal edge-on 
(top-left and bottom-right) projections. 

star formation are turned on the stellar disc is highly warped 
but by 2.5 Gyr it settles into a single plane. Figure [9] shows 
that by 6.5 Gyr (where the simulation ends) a rapidly ro- 
tating thin stellar disc supporting spirals has formed. Gas 
cools onto the stellar disc with a small (~ 10°) warp. Figure 
1101 shows the evolution of the disc orientation. The stellar 
disc never settles into an intermediate-axis orientation; at 
3 Gyr the disc is inclined by ~ 30° to this axis, increasing 
to ~ 55° by 6 Gyr. Thus even with the global gas angular 
momentum aligned with the intermediate axis, the disc can- 
not form in an intermediate-axis orientation even though the 
halo is only very mildly triaxial, with T ^ 0.93 throughout 
the inner 100 kpc before the disc forms. 



4 DISCUSSION 

We have shown that a disc can never remain with its mi- 
nor axis aligned with the intermediate axis of a triaxial halo 
(an "intermediate-axis orientation"). This is shown in a dif- 
ferent way in Figure 1111 which plots the evolution of the 
angle between the stellar disc angular momentum and the 
halo's intermediate axis. In models lAl, IA2 and LCI this 
angle increases rapidly until the disc is nearly orthogonal. 
In model Gil the disc is initially chaotic, but once it settles 
after 2.5 Gyr the angle increases throughout. This happens 
even if the disc cocoons itself by flattening the inner halo 
such that the minor axis of the net potential is perpendic- 
ular to the disc where it resides. Such a vertically flattened 
inn er halo is the expected configuration within 20 kpc for 
the lLaw fc Maiewskil (|2010l ) triaxial model of the Milky Way 
(Johnston private communication). In that case, the orbits 
of stars in the disc are stable. However, a disc grown in 
an intermediate-axis orientation gives rise to an instability 
in the halo. As a result the inner halo tilts rapidly (within 




Figure 10. Briggs figure for run Gil at 0.5 Gyr intervals. Dot- 
ted circles are spaced at 20° intervals, with the outer solid circle 
corresponding to 8 = 100° . The centre of the disc is indicated by 
the (red) open circle while the remaining disc annuli are indicated 
by (red) crosses. The open (blue) star, square and triangle sym- 
bols indicate the direction of the pre-disc halo long, intermediate 
and short axes, respectively. 




2 4 6 8 10 



t [Gyr] 

Figure 11. Tilting of the models away from the intermediate 
axis orientation. Different simulations are shown by different line 
styles as indicated. During the first 2 Gyr model Gil is highly 
warped before it settles into a coherent plane. 
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~ 4 Gyr), shepherding the disc along with it. A hallmark of 
this instability is that the disc tilts without precessing, as it 
stays near equilibrium with the tilting inner halo. 

We also showed, by means of a simulation with gas and 
star formation, that even if the gas angular momentum is 
along the intermediate axis, that the disc which forms is 
not in the intermediate-axis orientation. This happens even 
when the halo is only very mildly triaxial: in model Gil the 
halo density has T ~ 0.93. We conclude that discs cannot 
form in an intermediate-axis orientation, and even if they 
were perturbed into such an orientation, they would not last 
long in it. 

The shape of the LCI potential varies strongly inside 
~ 15 kpc but this part of the potential is poorly con- 
strained by the Sagit tarius Stream (but see, for instance, 
iLoebman et al.l I2OI2I . for other constraints). Beyond this 
radius, the potential shape varies quite slowly. The ratio 
(6/a)$ > 0.9 which is not much different from the LMIO 
model, wh ile (c/a)<s> ~ . 9, wh ich is larger than in the 
LMIO and iDeg fc Widrowl (|2013l 'l models. Thus the poten- 
ti al in model LCI i s clos er to spherical than the m odels 
of iLaw fc Maiewskil (|2010l ) and Dog & Widrovj (|2013l V The 
top and bottom panels of Figure [7] also shows the potential 
shape in models lAl and IA2. Both of these are quite pro- 
late, with model lAl closer to spherical than the best Milky 
Way model of DW13. These less spherical Milky Way mod- 
els would therefore probably also be highly unstable. Since 
the Milky Way has not experienced large interactions in the 
past few gigayears that might have put it in an intermediate 
axis orientation, it has had ample time to tilt, so it is very 
unlikely to be in such an orientation. 

Further difficulties for the Sagittarius Stream mod- 
els come from their failure to match the leading arm of 
the Stream well, and to produce the observed bifurcation 
(Bclokurov ct al. 2006). We note that the best-fit model of 
[Law fc Maiewskil 12010), while it does an excellent job of 
fitting much of the observational data, still has = 3.4 
(but in comparison, their spherical halo has — 9). In 
the past solutions of these problems have been sought in 
detai ls of the Sagittarius dwarf itself (e.g. Pcfiarrubia ct al. 
I2OIOI ). Here we have shown that triaxial models of the Milky 
Way which consistently find the disc in an intermediate axis 
orientation themselves can be ruled out. What then is the 
most promising way to improve Milky Way halo models of 
the Sagittarius Stream? 

The assumption of a constant shape within the region 
of the Sagittarius stream is unlikely to be correct; however 
halo shapes generally change sufflciently slowly beyond the 
disc that this assumption amounts to measuring an average 
shape rather than completely invalidating past models (note, 
for instance, how small the variation in the shape of the 
pot ential of model LCI is from 16-60 kpc in Figure [71 . 

Ilbata et al.l l|2012t ) showed that if the halo rotation 
curve is allowed to increase to ~ 300 kms~^ at 60 kpc that 
it is still possible to fit the Sagittarius stream by a spherical 
model. This model still fails to produce a bifurcation and 
resuhs in a quite massive Milky Way (2.6 - 3.1 x 10^^ M©). 
As argued byjbata et al. (20121), such a model cannot be ex- 
cluded by current observational constraints but it would be 
unusual in ACDM. Nonetheless more general density profiles 
are certainly highly recommended for future models. 

We propose here a different, and more natural, solution 



to the problem s of th e Sagittarius Strea m. The models o f 
iDeg fc Widr'wl l|2013h as well as those of iLaw et all (|2009l ) 
vary the axes ratios of the halo such that if the disc had been 
perpendicular to either the short or the long axes of the halo 
then the models would have been able to recover this; the 
fact that they did not means that the Milky Way disc is 
not in either orientation. We contend that the assumption 
that the disc of the Milky Way is in one of the symmetry 
planes of the halo must be incorrect. The possibility that this 
assumption can fail is clearly illustrated by our model Gil 
which shows that the disc does not need to be sitting in one 
of the principal planes of a triaxial halo outside the region 
dominated by the disc. Indeed in cosmological simulations a 
decoupling between t he disc /inner halo and the outer halo is 
a common outcome (|Bailin et al. I I2OO5I : iRoskar et al.ll2O10l '). 
Future models of the Sagittarius Stream need to consider 
the possibility that the disc is not in one of the symmetry 
planes of the halo. 

While complicating efforts at understanding the halo, 
this orientation nonetheless provides a unique op portunity 
to test the Modified Newtonian Dynamics (MOND: lMilgroml 
[l983; Bckciistein 2004) . If the Sagittarius Stream requires a 
net potential that is tilted with respect to the Milky Way 
disc, as we have argued, then this would constitute a prob- 
lem for MOND, which requires the symmetry axis of the 
disc and of the net potential to be parallel. The forthcom- 
ing ge neration of Milky Wa y surveys and missions such as 
Gaia (iPerryman et al.ll200ll) and the Large Synoptic Survey 
Telescope ( Ivezic et al.l " 20081 ') will provide the data needed 
for much more accurate modelling of the Milky Way's po- 
tential. 
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APPENDIX A: AN INTERPRETATION OF THE 
HALO INSTABILITY 

Here we explore the cause of the halo instability which pre- 
vents discs from inhabiting an intermediate axis orienta- 
tion. As shown above, the orientation of the inner potential 
changes as the disc is grown within it. In model lAl, the 
axes of the potential are initially ordered as > > J/*, 
but once the disc grows, the inner potential gets flattened 
and has axes ordered as x<s, > y<s, > z^, while at larg e r radi i 
the original axes ordering is retained. IValluri et al. (|201Cll ) 
showed that while tube orbits are uncommon in halo A be- 
fore the disc forms, a fraction of halo box orbits are trans- 
formed by the growing disc, with short axis tubes becom- 
ing abundant (we refer to the axes ordering at large radii, 
rather than in the flattened inner halo, to define orbit fam- 
ilies). Because of the radial change in the axes ordering, 
particles circulating about the short axis of the inner halo 
are actually intermediate axis tubes (lATs) if they venture 
outside the inner halo. We propose that tube orbits crossing 
the inner halo are destabilised by the radially varying halo 
orientation and drive the instability of the inner halo. We 
explore this hypothesis by comparing models lAl and LAI . 
Model LAl, which was presented by IValluri et al.l (|2012l ). 
is identical to model lAl other than that the disc is grown 
perpendicular to the long axis, which we found is a stable 
orientation for this disc. In model LAl, the original axes or- 
dering is 2:$ > > y<s becoming, in the inner (< 15 kpc) 
halo, > j/4> > z,^ once the disc is grown. As with lAl, 
any particles on tube orbits can be destabilised by crossing 
from the flattened inner halo to larger radii. Thus model 
LAl acts as a control in the interpretation of why lAl (and 
the intermediate-axis orientation in general) is unstable. 

Orbits of dark matter particles in LAl and lAl 
were an a lysed using the Laskar fr e quenc y analysis method 
iLaskail Il993l : IValluri fc MerrittI Il998h with the auto- 
mated orbit class i ficatio n scheme described previously 
IValluri et al.ll2010l . [20ll 'l. Briefly, Laskar's method uses a 
filtered Fourier transform method to obtain accurate orbital 
frequency spectra from complex time series constructed from 
the orbital phase space coordinates. The frequency spectra 
are then decomposed into the set of three linearly inde- 
pendent base frequencies (the "fundamental frequencies") 
of which all other frequencies in the spectrum are integer 
linear multiples. The ratios of fundamental frequencies are 
rationalised following a meth od similar to that described by 
Carpintero fc Aguilaj (| 19981 ). Previously (e.g. IValluri et al.] 



Model 


Total 


Boxes 


LATs 


SATs 


lATs 


lAl 


6697 


4157 


1400 


378 


762 


LAl 


6782 


3443 


1316 


2023 






2010l ) we only considered classification into the traditional 



orbit families believed to constitute triaxial galaxies (boxes, 
long-axis tubes, short-axis tubes, and various families of res- 
onant orbits). Here we adapted our code to also consider the 
possibility that orbits may be tubes which circulate about 
the intermediate axis. 

We measure the degree of diffusivity of individual or- 
bits via the diffusion rate parameter log(A/). Since regular 
orbits have fixed frequencies, a chaotic orbit can be identi- 
fied if its fundamental frequencies measured in the two con- 
secutive time segm ents change significantly (Laskar 199^. 
IValluri et al] (|2010l ) showed that even for orbits in A^-body 
potentials (which are inherently noisy) it is possible to dis- 
tinguish between A'^-body jitter and true chaos via a quan- 
titative measurement of frequency drift by defining log(A/) 



Table Al. The number of orbits in the different families in the 
two models from a sample of lOK orbits. Only those orbits which 
complete 30 periods in 50 Gyr integrations are counted. 



as the logarithm of the change in the frequency of the lead- 
ing term in the or bit's frequency spectrum in two consecu- 
tive time segments. I Valluri et al.l ( |2010 ) showed, using orbits 
in A'-body simulations of spherical halos, that orbits with 
log(A/) < —1.2 were regular. 

W e use the orbit sample described in IValluri et al.l 
(2012'): briefly, this is a sample of orbits for 10* particles 
in each model. Each of these particles was chosen at ran- 
dom from those within 200 kpc from the centre before the 
disc was grown; the same set of particles are used in mod- 
els lAl and LAl. Each orbit is integrated for 50 Gyr. The 
frequency analysis is not guaranteed to produce accurate fre- 
quencies for orbital integration times less than 20-30 orbital 
periods. Table lXTI lists the number of orbits of different types 
with more than 30 orbital periods in our sample. About two- 
thirds of all orbits satisfy the orbital periods condition; more 
than half of these are box orbits. Model lAl contains ~ 30% 
more box orbits than model LAl. This probably contributes 
to making it more unstable since box orbit have zero average 
angular momentum making them easier to tilt. 

Figure lAll plots the distribution of log(A/) for tube 
orbits of all types in models lAl and LAl, separated into 
three groups by radial 

5 kpc < rperi < rt < Vapo, and rt < r^eri, where r^eri 
and Tapo are the peri- and apo-centre distances and rt is 
the radius at which the potential switches orientation. From 
Figure|4]we find rt = 25 kpc for model lAl, whereas a sim- 
ilar measurement for LAl gives rt — 15 kpc. Orbits that 
never visit the inner region have low log(A/). Orbits that 
remain wholly within the inner region have higher diffusion 
rates, but they tend to be less numerous. Orbits that move 
across rt are the most abundant and have higher diffusivity 
in model lAl than in LAl. It is this difference in the diffu- 
sion of tube orbits crossing the radius at which the potential 
reorients that we propose accounts for the different stability 
properties of models lAl and LAl. 
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Figure Al. The distributions of log(A/) for tube orbits in model 
lAl (black lines) and LAI (red lines). The solid, dashed and dot- 
ted lines show those orbits with 5 kpc < Tpg,.i < i"apo < ft , 
5 kpc < Tpgri < rt < Tapo and Vp^^i > 25 kpc, respectively. For 
model lAl, rt = 25 kpc while for LAI rt = 15 kpc. 



